Test paleoTS with Fossil Checklist data (but is probably of no use, because they report average body sizes (means, median, something else? what are the respective sample size? maybe ask the authors!?), so this is just for playing around).
Raw data:
library(paleoTS)
#setwd("//naturkundemuseum-berlin.de/MuseumDFSRoot/Benutzer/Julia.Joos/Eigene Dateien/MA")
test<-read.csv("test26.5.csv", sep=";", header=TRUE)
test
The first plot shows mean Cl size for each taxon as a single data point, so each data point is one species (in this case this equals one individual, since I don’t have sample sizes), even within time bins.
Test1 <- test %>%
mutate(mm = CL_mean, vv=0, nn= n, tt=Age_mean) %>%
dplyr::select(mm, vv, nn, tt)
paleoTest1 <-as.paleoTS(Test1$mm, Test1$vv, Test1$nn, Test1$tt, MM = NULL,
genpars = NULL, label = "Testudinidae body size evolution mode")
paleoTest1
$mm
[1] 107.5 24.0 28.0 23.0 165.0 80.0 120.0 175.0 11.0 90.0 33.0 150.0 115.0 125.0 47.0 75.0 37.5 192.5 195.0 9.0
[21] 83.0 95.0 40.0 58.0 34.0 29.0 20.0 23.0 45.0 60.0 50.0 100.0 190.0 22.0 100.0 90.0 195.0 120.0 20.0 12.0
[41] 22.5 15.0 186.0 62.5 87.5 85.0 50.0 20.0 25.0 13.0 27.5 52.0 122.5 80.0 37.5 150.0 79.0 98.0 96.0 88.0
[61] 90.0 100.0 46.0 100.0 120.0 110.0 40.0 150.0 90.0 120.0 27.5
$vv
[1] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
[63] 0 0 0 0 0 0 0 0 0
$nn
[1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
[63] 1 1 1 1 1 1 1 1 1
$tt
[1] 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 -1.28765
[14] -1.28765 -1.28765 -1.28765 -1.28765 4.99550 3.23250 1.40950 -1.22465 -1.22465 -1.22465 -1.22465 -1.22465 -1.22465
[27] -1.22465 -1.23000 -1.23000 -1.23000 -1.23050 0.56950 0.56950 0.56950 0.56950 0.56950 0.56950 0.56950 0.56950
[40] 0.56950 0.56950 0.56950 0.56950 0.56950 -0.84000 -0.84000 -0.84000 -0.84000 -0.84000 -0.84000 -0.84000 -0.84000
[53] -0.89715 -0.89715 -0.89715 -0.89715 -1.29350 -1.29350 -1.29350 -1.29350 -1.29350 -1.29350 -1.29350 -1.29350 -1.29350
[66] -1.29350 0.00635 0.00635 0.05050 1.00650 1.00650
$MM
NULL
$genpars
NULL
$label
[1] "Testudinidae body size evolution mode"
$start.age
[1] 1.2935
$timeDir
[1] "increasing"
attr(,"class")
[1] "paleoTS"
plot(paleoTest1)
This is the underlying data for Test1:
Test1
For the second plot, I averaged CL means across taxa for each time bin, which leaves one data point per time bin, comprising all taxa within the respective bin:
Test2 <- test %>%
group_by(Age_mean) %>%
summarise(mm = mean(CL_mean), nn=n(), vv=var(CL_mean)) %>%
mutate(tt=Age_mean) %>%
dplyr::select(mm, vv, nn, tt)
# NA: column 2, rows 3, 10, 13, 14, 15
Test2[3,2] <- 0
Test2[10,2] <- 0
Test2[13,2] <- 0
Test2[14,2] <- 0
Test2[15,2] <- 0
paleoTest2 <-as.paleoTS(Test2$mm, Test2$vv, Test2$nn, Test2$tt, MM = NULL,
genpars = NULL, label = "Testudinidae body size evolution mode")
paleoTest2
$mm
[1] 92.70000 79.90000 50.00000 42.66667 51.28571 97.50000 45.00000 83.87500 95.00000 90.00000 87.30769 73.75000
[13] 9.00000 195.00000 192.50000
$vv
[1] 398.6778 1542.5500 0.0000 346.3333 810.5714 2429.1667 833.6429 3589.5511 6050.0000 0.0000 4816.0224 4278.1250
[13] 0.0000 0.0000 0.0000
$nn
[1] 10 5 1 3 7 4 8 12 2 1 13 2 1 1 1
$tt
[1] 0.00000 0.00585 0.06300 0.06350 0.06885 0.39635 0.45350 1.29350 1.29985 1.34400 1.86300 2.30000 2.70300 4.52600 6.28900
$MM
NULL
$genpars
NULL
$label
[1] "Testudinidae body size evolution mode"
$start.age
NULL
$timeDir
[1] "increasing"
attr(,"class")
[1] "paleoTS"
plot(paleoTest2)
Since “real” variances and sample sizes are available when pooling all taxa, you can even fit models (as you should be able to in the end). (when I remember correctly, the model with the highest Akaike.wt is the best supported one, in this case this would be URW = random walk)
a=fit3models(paleoTest2, silent=FALSE, method="AD", pool=FALSE) #not working with Test1, because no variances/sample sizes available, I guess
Comparing 3 models [n = 14, method = AD]
logL K AICc Akaike.wt
GRW -70.40398 2 145.8989 0.373
URW -71.26818 1 144.8697 0.625
Stasis -75.70460 2 156.5001 0.002
str(a)
'data.frame': 3 obs. of 4 variables:
$ logL : num -70.4 -71.3 -75.7
$ K : num 2 1 2
$ AICc : num 146 145 157
$ Akaike.wt: num 0.373 0.625 0.002
a$AICc[1] # not sure what this tells me...
[1] 145.8989
This is the underlying data for Test2:
Test2
Try paleoTS with some first real data. Here is the underlying data:
tidyCL<-read.csv("tortoises_tidy.csv", sep=";", header=TRUE)
tidyCL
Prepare data for conversion to paleoTS-object:
TidyCL <- tidyCL %>%
select(MAmin, Mamax, CL) %>%
filter(CL != "NA") %>%
mutate(tt= (MAmin+Mamax)/2) %>% # create mean age
group_by(tt) %>% #create time bins
summarise(mm=mean(CL), vv=var(CL), nn=n()) #create means etc. for each time bin
TidyCL[is.na(TidyCL)]<-0 #subset NAs with O for
TidyCL
bins <- tidyCL %>%
# select(MAmin, Mamax, CL) %>%
filter(CL != "NA") %>%
mutate(tt= (MAmin+Mamax)/2) %>% # create mean age
group_by(tt)
bins
library(paleoTS)
paleoTidyCL <-as.paleoTS(TidyCL$mm, TidyCL$vv, TidyCL$nn, TidyCL$tt, MM = NULL, genpars = NULL, label = "Testudinidae body size evolution mode")
paleoTidyCL
$mm
[1] 195.0000 180.5000 860.3333 880.0000 1200.0000 278.0000 1200.0000 500.0000 500.0000 190.7500 370.0000 400.0000
[13] 400.0000
$vv
[1] 0.0000 40.5000 307760.3333 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 911.9286
[11] 0.0000 0.0000 0.0000
$nn
[1] 1 2 3 1 1 1 1 2 1 8 1 1 1
$tt
[1] 0.000 1.230 2.180 2.730 3.730 6.730 7.730 8.330 10.080 11.230 16.730 21.345 32.180
$MM
NULL
$genpars
NULL
$label
[1] "Testudinidae body size evolution mode"
$start.age
[1] 1.77
$timeDir
[1] "increasing"
attr(,"class")
[1] "paleoTS"
plot(paleoTidyCL)
fit3models(paleoTidyCL, silent=FALSE, method="AD", pool=FALSE) #not working with Test1, because no variances/sample sizes available, I guess
Comparing 3 models [n = 12, method = AD]
logL K AICc Akaike.wt
GRW -94.17833 2 193.6900 0.001
URW -104.38851 1 211.1770 0.000
Stasis -87.43929 2 180.2119 0.999
Map <- tidyCL %>%
select(Genus, Taxon, Latitude, Longitude, Country, CL, PL) %>%
group_by(Latitude) %>%
mutate(count= n())
mapWorld <- borders("world", colour="azure3", fill="azure3") # create a layer of borders
mp <- Map %>%
ggplot(aes(Longitude, Latitude)) + mapWorld +
# geom_point(fill="red", colour="red", size=0.5) +
geom_point(aes(Longitude, Latitude,colour=CL, size=count))
mp
library(plotly)
ggplotly(mp)
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
Map all localities with sample size and age indicated (regardless of whether CL information is available):
test<-read.csv("tortoises13-04.csv", sep=";", header=TRUE)
colnames(test)[6] <- "Mamin"
colnames(test)[7] <- "Mamax"
Test <- test %>%
select(Locality, Country, Latitude, Longitude, Mamin, Mamax, Epoch, Genus, Species, Taxon, CL) %>%
mutate(Age= (Mamin+Mamax)/2) %>% # create mean age
group_by(Latitude) %>%
mutate(count= n())
mapWorld <- borders("world", colour="azure3", fill="azure3") # create a layer of borders
map <- Test %>%
ggplot(aes(Longitude, Latitude)) + mapWorld +
#geom_point(fill="red", colour="red", size=0.5) +
geom_point(aes(Longitude, Latitude,colour=Age, size=count))
map
ggplotly(map)
We recommend that you use the dev version of ggplot2 with `ggplotly()`
Install it with: `devtools::install_github('hadley/ggplot2')`
This is an R Markdown Notebook. When you execute code within the notebook, the results appear beneath the code.
Try executing this chunk by clicking the Run button within the chunk or by placing your cursor inside it and pressing Ctrl+Shift+Enter.
Add a new chunk by clicking the Insert Chunk button on the toolbar or by pressing Ctrl+Alt+I.
When you save the notebook, an HTML file containing the code and output will be saved alongside it (click the Preview button or press Ctrl+Shift+K to preview the HTML file).